Universality class of the pair contact process with diffusion 
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The pair contact process with diffusion (PCPD) is studied with a standard Monte Carlo approach 
and with simulations at fixed densities. A standard analysis of the simulation results, based on the 
particle densities or on the pair densities, yields inconsistent estimates for the critical exponents. 
However, if a well-chosen linear combination of the particle and pair densities is used, leading 
corrections can be suppressed, and consistent estimates for the independent critical exponents 5 = 
0.16(2), [3 — 0.28(2) and z = 1.58 are obtained. Since these estimates are also consistent with their 
values in directed percolation (DP), we conclude that PCPD falls in the same universality class as 
DP. 

PACS numbers: 05.70.Ln, 05.50.+q, 64.60.Ht 



I. INTRODUCTION 

In the fermionic one-dimensional pair contact process 
with diffusion (PCPD), the model studied here, point 
particles can diffuse on a line, and reactions can occur 
when two particles end up next to each other. The par- 
ticles can then both be annihilated, or if there is a free 
site next to the pair, a new particle can be created. The 
reactions present in the PCPD model are: 
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with < d < 1 and < p < 1. 

As for other systems of this kind, the scaling relations 
which are expected to hold for the density p and correla- 
tion length £ are 
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in which e = p c — p is the distance to the critical point. 
The critical exponents in these scaling relations define 
the universality class a system belongs to[l|, [2|. Two 
firmly established universality classes so far for systems 
of this kind are the Directed Percolation (DP) and Parity 
Conserving (PC) classes. 

Since the introduction of PCPD a decade ago or at 
least a model closely resembling it, it has attracted much 
attention. The main reason for this attention is that, 
while its symmetries and conservation laws are seem- 
ingly identical to directed percolation (DP), numerical 
estimates of its critical exponents seem to place it in a 
different universality class. It has been conjectured be- 
fore by Grassberger |i[ and Janssen 0] that critical points 
with a unique absorbing state and a single order param- 
eter will fall into the DP universality class, making this a 



likely candidate for PCPD as well. However, the fact that 
the PCPD absorbing state is not unique makes this less 
certain, and it has been suggested by earlier studies that 
it might not be possible to describe PCPD with a single 
order parameter [HQ- Table Q] shows values reported for 
the critical exponents of PCPD from previous studies, as 
well as the accurately known values for these exponents 
in the DP universality class. So far, while it is obvious 
much research has been done on this subject, there is still 
a certain amount of disagreement in the results. Gener- 
ally, the differences between the measured exponents and 
the values for directed percolation have been significant. 
It has been argued, however, that the discrepancy in S 
between PCPD and DP is due to severe finite-size and 
finite-time effects Q , and recent ultralong simulations @ 
show a clear trend of 5 towards its DP-value. 

Other recent studies have looked at field theory ap- 
proaches for PCPD Q , and looked at a bosonic version of 
the model, where multiple particles can exist in one spot 
0, [1(3] ■ Another study (Tl| has investigated the structure 
and behavior of clusters in PCPD, to clarify the slow 
approach of PCPD to its asymptotic scaling regime. In- 
vestigation of the crossover from PCPD to DP [l2| has 
yielded evidence towards non-DP scaling. 

At this point, it is still unclear what universality class 
PCPD does belong to. Most studies so far conclude that 
PCPD might belong to a new univerality class [l(| [12], 
or even that it might belong to several ones based on 
the value of the parameter d fl3l [Tl . flU fTsj ] . On the 
other hand, in a recent study by H. Hinrichsen the 
critical exponent 5 was shown to display a significant 
drift towards the DP value, providing evidence towards 
a single universality class. However, this tendency has 
not been clearly shown for the other two exponents so 
far. 

A previous study (Ref. 0]) has provided numerical 
evidence that, at the critical point, the ratio between 
the pair and particle densities tends to a non-zero, finite 
value when the simulation time tends to infinity. Also the 
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high-quality data of Ref. @ confirms the convergence to 
a nonzero ratio. A visual inspection of the system shows 
big clusters, separated by increasingly large empty (or 
near-empty) regions. Since these clusters tend to have a 
finite density in particles as well as in pairs, the conver- 
gence of their ratio to a nonzero value is not unexpected. 

Based on this convergence of the ratio of particles and 
pairs, we will show in this paper that PCPD has strong 
correction terms in its scaling relations when it has not 
yet reached the thermodynamical limit. These correction 
terms will distort any direct measurements of the expo- 
nent. The correction terms for these quantities, however, 
do not turn out to be exactly the same. We will show 
that if both the particle- and the pair-densities are com- 
bined in the analysis, these corrections can in some cases 
be suppressed, allowing a more accurate result for the 
exponents. The resulting estimates for the critical ex- 
ponents of the PCPD model are consistent with the DP 
universality class. 



II. SIMULATION APPROACH 

In our analysis, we use both normal Monte Carlo sim- 
ulations and simulations at constant density. For the 
former, we followed the usual method of simulating a 
system of sites that can each contain a particle, using 
the same multispin coding program as used by Barkema 
and Carlon in an earlier study @. 

It is very hard to get accurate data at low densities 
(below 10%) with standard simulations, due to time con- 
straints and the possibility of the system reaching an ab- 
sorbing state because of its finite size. Thus, to get more 
accurate data at lower densities, we also performed sim- 
ulations of the PCPD at constant density. This was done 
by having two possible reactions for the system: each step 
consists of cither the usual diffusion, or a combination of 
one annihilation of a pair of particles and the creation 
of two new ones. It has been shown [23|, [U HH that 
this procedure, after thermalization, produces configu- 
rations which are indistinguishable from those obtained 
with standard simulations (with a fixed value for p < p c ) 
after the same density has been reached. 

To correctly update the rates of the reactions in this 
system, the number of possible places where these reac- 
tions can take place is required. Counting these at each 
time step would take too much time, so instead of keeping 
track of sites, and whether or not they are occupied, only 
the existing particles and the distance to their neighbors 
are stored during the simulation. Based on the number of 
direct neighbors of a particle (0, 1 or 2), the information 
is kept in one of three lists. Particles are moved from one 
list to another when a reaction causes them to gain or lose 
a neighbor. This method allows us to track the number 
of particles with 0,1, and 2 neighbors at every step of the 
simulation. In addition, it enables us to pick randomly a 
particle with a certain number of neighbors, making sure 
that we can always perform a chosen reaction. We can 



just compute the probabilities for each reaction in a given 
configuration, choose a reaction based on that, and pick 
a site to perform it. This prevents the rejected reactions 
normally seen in Monte Carlo simulations, allowing us to 
speed up the simulations significantly. 

The probability to perform a reaction in a general 
PCPD system simply equals the rate (same as in eq. (fT])) 
multiplied by the number of possible places where the 
reaction can be performed, yielding 
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Here, rij is the number of particles with i direct neigh- 
bours, and the probability for each reaction is not nor- 
malized yet. Since we want the probability for procre- 
ation to be twice as large as that for annihilation to keep 
a constant number of particles, we can compute the value 
of p that would achieve this from the values of n^ at each 
simulation step. Knowing this value at each step allows 
us to compute the probabilities for each reaction in the 
model with constant density using 
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Only two reactions exist now: diffusion, and the combi- 
nation of two particle creations and one pair annihilation. 
Since the latter reaction is actually a combination of three 
separate events, the probability for diffusion is multiplied 
by a factor of three to compensate this. The computed 
value for p e g can be monitored over the course of the sim- 
ulation and averaged to find the value of p corresponding 
to the fixed density of the system. 

Of course, quantities depending on time cannot be 
measured from these simulations at constant density. 
However, these simulations are very useful for determin- 
ing the exponent /?, since it allows us to explore the re- 
lation between p and p for low densities much faster and 
without the risk that the system reaches an absorbing 
state by fluctuations. 

Apart from the usual density, the pair density p* was 
also measured in both normal and constant-density sim- 
ulations. This density is simply the number of pairs of 
directly adjacent particles, divided by the length of the 
system. Since the ratio between p and p* approaches a 
non-zero constant in the thermodynamic limit, as shown 
by Barkema and Carlon [8], p* should obey the same 
power laws as p, with the same exponents. Finite-size or 
finite-time effects might however not be exactly the same 
for p and p* . Therefore, if we find different critical expo- 
nents for the particle-density and pair-density, we know 
that the method of analysis used is incorrect. 
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Study 


Year 


d 


5 





z 




Odor [16] 


2000 


0.1 
0.5 
0.9 


0.275(4) 

0.21(1) 

0.20(1) 


0.58(1) 
0.40(2) 
0.39(2) 


- 


- 


Carlon, Henkel 
and Schollwock [17] 


2001 


0.1 
0.5 
0.8 


- 
- 


- 
- 


1.87(3) 
1.70(3) 
1.60(5) 


0.50(3) 
0.48(3) 
0.51(3) 


Hinrichsen [181 


2001 


0.1 


0.25 


< 0.67 


1.83(5) 


0.50(3) 


Park, Hinrichsen & Kim [19] 


2001 




0.236(10) 


0.50(5) 


1.80(2) 




Park & Kim [13] 


2002 




0.241(5) 
0.242(5) 


0.496(22) 
0.519(24) 


1.80(10) 
1.78(5) 




Dickman & de Menezes [14] 


2002 


0.1 
0.5 
0.85 


0.249(5) 
0.236(3) 
0.234(5) 


0.546(6) 
0.468(2) 
0.454(2) 


2.04(4) 
1.86(2) 
1.77(2) 


0.503(6) 
0.430(2) 
0.412(2) 


Odor [15] 


2003 


0.1 
0.5 
0.7 


0.206(7) 
0.206(7) 
0.214(5) 


0.407(7) 
0.402(8) 
0.39(1) 


1.95(1) 
1.84(1) 
1.75(1) 


0.49(2) 
0.41(2) 
0.38(2) 


Kockelkoren & ChatcflO] 


2003 


* 


0.200(5) 


0.37(2) 


1.70(5) 




Barkema & Carlon [8] 


2003 


0.1 
0.2 
0.5 
0.9 


0.17 
0.17 
0.17(1) 
0.17 


- 
- 
- 


- 

1.70(1) 
- 

1.61(3) 


- 

0.28(4) 
0.27(4) 


Noh & Park [20] 


2004 


0.1 


0.27(4) 


0.65(12) 


1.8(2) 


0.50(5) 


Park & Park [7] 


2005 


1/3 


0.20(1) 








Hinrichsen [9] 


2006 


* 


< 0.185 


< 0.34 


< 1.65 


- 


De Oliveira & Dickman [21] 


2006 


0.1 
0.5 
0.85 


- 


- 


2.08(15) 

2.04(5) 

1.88(12) 


0.505(10) 
0.385(11) 
0.386(5) 


Kwon & Kim [11] 


2007 








1.61(1) 




Directed Percolation 






0.1595 


0.2765 


1.5807 


0.2521 



TABLE I: Reported values for the critical exponents of PCPD [22] . In some of the studies, a modified model was used, 
changing the definition of the parameters; these studies are marked with a star in the d column. 



III. DIRECT COMPUTATION OF CRITICAL 
EXPONENTS 

A direct analysis of the simulation data from all per- 
formed simulations is fairly straightforward, and has al- 
ready been shown before to give rise to great inaccuracies 
in its results. In this section we will give a quick overview 
of this direct analysis, with particular attention for the 
differences between the results for the particle density 
p and pair density p*. A first estimate of the critical 
exponent 8 is obtained by taking the logarithm of both 
the density and time in a simulation close to the criti- 
cal point, and fitting a straight line through the data. 
The simulations used to determine 8 were run on a sys- 
tem with L — 100 000, for about 3 • 10 6 timesteps. For 
d = 0.5, this leads to 8 = 0.19, as shown in Figure [T] A 
slightly higher value for 5 is found, however, if this anal- 
ysis is performed on the pair density instead (8* = 0.20), 
and it is visible that the double-logarithmic curve is not 



entirely straight. As was already shown before d,[{|, this 
curving tendency can very well be extrapolated to <5dp- 

The value of f3 can be determined by performing sim- 
ulations at non-critical values of p, where the system 
reaches its steady state. Here, we used simulations at 
constant density for this. The exponent j3 can be ex- 
tracted from a logarithmic plot of the density versus 
p — p c . Figure [5] shows such a plot. At higher densi- 
ties, this yields values for /3 which are consistent with 
earlier reported measurements, significantly higher than 
/3dp- However, it is clear that this fit does not hold 
up for the lower densities, which can be reached in the 
constant-density simulations. In addition, for the region 
with higher densities, the slopes for the particle and the 
pair density are clearly different, showing that an analysis 
consisting of a simple linear fit cannot be trusted. 

For the calculation of the critical exponent z = vh/v± 
we examine finite-size effects on a system at the critical 
value of p. If the system size L is small, these effects will 
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FIG. 1: Particle density p and pair density p* as a function 
of time, for PCPD simulations at d = 0.5, p = 0.1524 and 
L — 100, 000, averaged over 64 simulations. The straight line 
is a fit to determine S ignoring correction terms. The top 
curve shows the density, with ^direct = 0.19. The bottom 
curve shows the pair density, with ^direct = 0.20. 
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FIG. 2: The particle density p and pair density p* as a func- 
tion of p — p c in constant-density simulations, with d = 0.5 
and L = 100, 000. The data for the pair density is the lower 
set. The fits, using only the points with higher densities, yield 
Aiircct = 0.42 for the particle density, and /3direct = 0.51 for 
the pair density. 

cause the density to start decaying exponentially once 
the correlation length £ ~ t 1 ^ approaches the system 
size. Since the particle and pair densities are tied closely 
together, they will collapse at the same time. After suf- 
ficient time, the density will then decay as: 

p ~ exp(-b(t/L z )) & (5) 

ln{p) = a-b{t/L z ). (6) 

Using Eq. (JBJ), the exponent z can be obtained from 
simulations in small systems of various sizes until past 
this point of collapse, and adjusting the exponent z un- 
til a data collapse can be attained for the exponential 
regime. We simulated systems for d — 0.1,0.5 and 0.9, 
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FIG. 3: Data collapse to determine z, ignoring correc- 
tion terms. Each curve is the average over 3200 simula- 
tions, for d = 0.5, p = 0.1524, and L = 200(+), 300 (x), 
500(*), 1000(H), 2000(H), 3000(0), and 5000 (•). we find 
Zdirect = 1.75. Since the particle density of the system will not 
always tend to due to the possibility of a single remaining 
particle, we only include the data on the pair density. 



d 


Pc 


6 


S* 
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f3* 


z 


0.1 


0.1111 


0.22 


0.24 


0.48 


0.60 


1.83 


0.5 


0.1524 


0.19 


0.20 


0.42 


0.51 


1.75 


0.9 


0.2333 


0.19 


0.20 


0.34 


0.39 


1.64 



TABLE II: The results of the direct analysis of the critical 
exponents, ignoring correction terms. 



with system sizes ranging from L = 200 to 5000 sites. 
For system sizes larger than this, the collapse occurs at 
late times, and therefore at a very low density. Apart 
from the very long simulation times, this poses another 
problem. To have accurate data at that point, we would 
have to know p c accurately enough such that its error 
Ap c obeys (Ap c )@ <C p. With the density dropping be- 
low 0.04 near the point of collapse for larger systems, Ap c 
would have to be smaller than 10~ 5 to avoid significant 
systematical errors. Since our precision in p c is not that 
accurate, our range for the system size is limited by this 
effect. 

At the time the density starts collapsing, t ~ L z . Since 
up to this point, the density was following a power law, 
we know that the density at this point will be p co u ~ 
t~ s ~ L~ Sz . Using this, we can scale the data from our 
simulations to obtain a data collapse, as shown in Figure 
[3] for d = 0.5. It turns out that scaling the vertical axis 
works best with Sz = S^pzop, though the optimal value 
for z in the horizontal scaling again varies with d. 

Table |TT] shows the results for all exponents, for the 
three values of d we investigated. The exponents all vary 
when the diffusion parameter d changes, which, assuming 
PCPD falls into a single universality class, again points 
out there is something wrong with such a direct analysis. 
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FIG. 4: The decay in time of linear combinations of the 
particle and pair densities, plotted to find the combination 
that isolates the correction term. Here, pe = a/a* p* — p, with 
a/a* = 2.41, 2.43 and 2.45 (from bottom to top), at d = 0.5 
and p = 0.1524. The top and bottom lines are shifted up and 
down by 1 for clarity. The middle plot is straightest; from its 
slope we get the correction exponent 9 = 0.63. 

IV. ANALYSIS INCLUDING CORRECTIONS 

Given the problems with the direct analysis of the crit- 
ical exponents 5 and j3, we propose introducing correc- 
tion terms into our scaling laws. Both p and p* obey 
the power laws as before, but with the DP values for the 
leading term in each relation, and an added correction 
term with a higher exponent. For the time-dependence 
of the density, this yields 

p = at~ 6 + bt~ e . , 

p* = a*t- s + b*t- e . 

While it is expected that the exponents for the two rela- 
tions are the same, the prefactors for the particle density 
can be different from those for the pair density. This 
makes it possible to determine linear combinations of p 
and p* where one of the two terms is suppressed, using 

p e = JLp*-p = ct °. 

To determine all coefficients accurately, we will need both 
exponents. The correction exponent 9 can be computed 
once a value for a/ a* has been found that turns hi pg into 
a linear function of lni, at least for low densities. As a 
starting estimate for this ratio, an extrapolation of the 
ratio p/p* for p — > can be used. Figure [5] shows the 
plots used to determine a/a* and calculate 9. 

With the exponent 9 known, and using the DP value 
for S, it is now possible to fit both p and p* with a lin- 
ear combination of t~ s and t~ e , to obtain the prefactors 
a, b, a* and b* . These fits are shown for d = 0.5 in Figure 
[5j The consistency of our fits can then be checked by 
determining 5 and 9 again from a linear fit to logarith- 
mic plots of the appropriate linear combinations of p and 
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FIG. 5: Time-dependance of the particle density p (top) 
and pair density p* (bottom) at d = 0.5, p = 0.1524 and 
L = 100, 000, averaged over 64 simulations. The lines were 
obtained from a least-squares fit of both t~ s and the correc- 
tion term t~ e to the data, using S = <5dp and 8 = 0.63. 



d 


a/a* 


9 


a 


6 


a* 


b* 


5 


Pcheck 


0.1 


2.56 


0.531 


0.41 


3.60 


0.16 


2.44 


0.16(2) 


0.571 


0.5 


2.43 


0.625 


0.47 


3.47 


0.19 


1.96 


0.16(2) 


0.633 


0.9 


4.19 


0.858 


0.37 


27.0 


0.089 


9.31 


0.16(2) 


0.887 



TABLE III: The results of our analysis for S, if a correction 
term is included. The # c hcck column shows the value of the 
exponent 9 obtained by using a and a* from the fit, to cross- 
check the values obtained in Figure [4] 



p*. In addition, the assumption for a /a* can be checked 
to make sure that it equals the ratio between the values 
from the fit. 

Going through this process for d — 0.5, the final fit 
to determine S from the computed ideal combination of 
p and p* gives a value of 0.16(2), as shown in Figure O 
consistent with Sbp = 0.1595. The Figure also shows the 
effective exponent as it changes during the simulation. 
The remaining curvature in these graphs is sensitive to 
small changes in the estimation of a/a*, even if those do 
not significantly affect the resulting exponents. There- 
fore, this deviation from a straight line is most likely 
caused by an inaccuracy in this estimation. The differ- 
ences between the estimated and calculated ratio a /a* 
and between the two obtained values for 9 are well within 
the error margins for those values. In Table IHIi the re- 
sults for the fit are shown also for other values of d. 

As seen in the table, the exponent 9 seems to vary as d 
changes. It is likely, however, that this is the effect of fur- 
ther correction exponents, whose coefficients depend on 
d. Given the numerical precision of our data, it would be 
too optimistic to claim that our values for 9 are accurate. 
A fit of Eq. (7) with 6 = 0.1595 to the high-quality data 
of Ref. @, which runs up to time t — 10 8 , yields correc- 
tion exponents as low as 9 = 0.3. This shows that either 
the corrections to 9 are very strong, or that the leading 
finite-time corrections actually cancel out in pg, causing 
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b 


a* 


b* 




Cchcck 


Pc 


0.1 


1.058 


0.62 


43.1 


0.24 


26.5 


0.29(3) 


1.055 


0.1111 


0.5 


1.056 


0.70 


16.4 


0.29 


10.5 


0.28(2) 


1.070 


0.1524 


0.9 


1.270 


0.43 


17.0 


0.10 


7.77 


0.29(3) 


1.285 


0.2333 



TABLE IV: The results of our analysis for (3, when including 
a correction term. The ("check column shows the value of the 
exponent £ obtained by using a and a* from the fit, which 
allows us to confirm our calculation of (. 
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FIG. 6: Fit to determine 8 from ps, the ideal linear com- 
bination of the particle and pair density that suppresses the 
correction term. The slope of the fit yields 8 = 0.16(1), at 
(from the top graph to the bottom) d — 0.1 and p = 0.1111, 
d = 0.5 and p = 0.1524, and d = 0.9 and p = 0.2333, with all 
data averaged over 64 simulations at each diffusion rate. The 
insets show the effective exponent <5 e ff = d\n(ps)/d\nt as a 
function of lnt, with the horizontal line at 8 — 0.16. 



us to measure the next correction exponent instead. 

A similar process can be followed for determining (3, 
using the constant density simulations. Again, we used 
L = 100, 000, with simulation times varying based on 
the relaxation time of the system. Our densities range 
from 0.05 to 0.4. For low densities, it can take up to 
10 9 simulation steps until the system no longer shows a 
systematic drift. The assumed behavior of the densities 
is 
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(9) 



Of course, the values for the prefactors, as well as the 
correction exponent, will be different from those in equa- 
tions [JJ However, the value for a/a* should still be the 
thermodynamic limit of the ratio p/p*, and thus will be 
the same as in our calculation of 6. 

Plotting the linear combination of p and p* which sup- 
presses the leading term, using the same ratio of a /a* as 
before, we can determine the first correction exponent £. 
With the exponents of the two leading terms known, we 
can again fit these to the data to determine their prefac- 
tors, and find the linear combination which will suppress 
the correction term. The result of this is shown in fig- 
ure [7J Again, the values from the fit to cross-check a /a* 
and £ show that these estimates are consistent with the 
fit. Tabic IIVI shows the results for different values of 
d. Again, the correction exponent varies slightly as the 
diffusion parameter d changes, suggesting additional cor- 
rections beyond the first. However, all data is consistent 
with a /3 DP = 0.2765. 

To confirm that we are performing our simulations at 
the critical point, we can calculate p c from these constant 
density simulations as well. With the same linear com- 
bination of the density and pair density as used for cal- 
culating (3, we again suppress the correction term. Since 
(pp) 1 ^ 3 is a linear function of p, and equals at p = p c , 
finding p c is a matter of linear extrapolation, as shown 
in the insets in Figure [JJ The results are in agreement 
with the values for p c used for all our simulations at the 
critical point, as seen in Table ITVl 

Lastly, we turn to the exponent z. While the data 
collapse in figure [3] is acceptable, we do find different 
values of Zdirect at d = 0.1,0.5 and 0.9. Since we have 
demonstrated that equal values for the exponents f3 and 
5 can be obtained for different values of d by including 
correction terms, this encourages us to try out a similar 
approach here. 

As before, we add a correction term to the relevant 
equation, in this case the one for the density during the 
collapse, and see how this fits with the DP exponents. 
The density will now decay differently, following: 



Pcoll 



exp(-&t(L^ + cL- x )), 



(10) 



with b an unknown constant and x t ne correction ex- 
ponent, with prefactor c. Considering that one of the 
determining factors of the system is diffusion, where the 
correlation length grows as y/t, a correction exponent of 
X = 2 is a reasonable assumption. Sadly, with the range 
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FIG. 7: Plots of the linear combination pp of the particle 
density and pair density that suppresses correction terms in 
the calculation of /3, using constant-density simulations at 
d — 0.1 (top), d — 0.5 (middle) and d — 0.9 (bottom), at 
L = 100, 000. The slope of the line gives (3 = 0.29(3), 0.28(2) 
and 0.29(3), from top to bottom, in agreement with DP. The 
insets show how p c for each value of d was checked using these 
simulations, by plotting the same {pp) 1 ^ as a function of p. 



of data available it is impossible to determine this expo- 
nent accurately. We can, however, assume z = zdp and 
X = 2, and show that this will provide data collapses at 
least as acceptable as the ones without a correction term. 
We use the same scaling as before, though we include the 
corrections calculated in our analysis for S in our vertical 
scaling: 



Pcoll = p*(t = L z ) 



■Sz 



b* 



(11) 



with 0, a* and b* taken from Table lllll 

Note that there still is only one free parameter in our 
data collapse after these assumptions have been made: 
instead of varying the exponent z, modifying the prefac- 
tor c is now the only way to make the data from different 
system sizes overlap, both for the particle density and the 
pair density. This prefactor for the correction term can 
vary as d changes, however. By plotting this relation for 
different sizes, we can find the value for c that causes a 
data collapse, as well as the equivalents for the pair den- 
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FIG. 8: Data collapses with z = zdp for d — 0.1, p — 0.1111 
(top), d = 0.5, p = 0.1524 (middle) and d = 0.9, p = 0.2333 
(bottom), at system sizes L = 200(+), 300(x), 500(*), 
1000(H), 2000(B), 3000(0) and 5000 («). All data is aver- 
aged over 3200 runs. The values for c used are c = 20, 15, and 
3, in order of increasing d. 



sity. The ratio of the particle density and pair density 
varies very little in the regime we are examining here, 
and the point of collapse is equal for both, so we expect 
that b and b* are equal. Fig. [5] shows the data collapses 
for three values of d. While this does show that all of 
our data can be seen as consistent with zdp, there is no 
clear way to show that this interpretation is better than 
the simpler approach, which yields varying exponents for 
different values of d. However, it does seem more likely 
that if the other two exponents are independent of the 
diffusion rate, the same should hold for z. 



V. CONCLUSION 

Monte Carlo simulations, both the usual approach with 
constant rate p and a new approach at constant density, 
have been used to analyze the critical behavior of the 
pair contact process with diffusion. While many recent 
studies conclude that this model does not belong to the 
directed percolation universality class, we find that if cor- 
rection terms are included in the power laws governing 
critical scaling, all of the acquired simulation data is con- 
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sistent with the exponents from DP. 

In our simulations, especially the calculations for (3 and 
5 offer convincing evidence that the DP values for these 
critical exponents are indeed accurate for the PCPD 
model. Since the critical exponents for p and p* must 
be equal, any linear combination of these must have the 
same exponents as well, in the thermodynamical limit, 
except in the singular case where the leading terms can- 
cel out. The fact that there exists a linear combination 
which follows a power law that is consistent with the DP 
exponent shows that the correct exponent for the system 



can at least not be any larger than that. Our analysis 
of z does not lend itself to an accurate calculation of the 
exponent, but does also show that a correction term can 
explain the deviation from the DP exponent. With all of 
our simulation data consistent with DP for all three of 
the studied exponents, and for all investigated diffusion 
rates, we conclude it is likely that PCPD does belong to 
the directed percolation universality class. 

We thank Haye Hinrichsen for providing us with the 
data described in Ref. @. 
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